home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Pascal Super Library
/
Pascal Super Library (CW International)(1997).bin
/
MATH
/
PMAT12
/
PMAT.EXE
/
PMATOP.PAS
< prev
next >
Wrap
Pascal/Delphi Source File
|
1993-01-17
|
48KB
|
1,555 lines
{
P-Mat - A Turbo Pascal program for Recursive Matrix Algebra
Version: 1.2
Author: Mark Von Tress, Ph.D.
Date: 01/30/93
Copyright(c) Mark Von Tress 1993
DISCLAIMER: THIS PROGRAM IS PROVIDED AS IS, WITHOUT ANY
WARRANTY, EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED
TO FITNESS FOR A PARTICULAR PURPOSE. THE AUTHOR DISCLAIMS
ALL LIABILITY FOR DIRECT OR CONSEQUENTIAL DAMAGES RESULTING
FROM USE OF THIS PROGRAM.
}
Unit pmatop;
Interface
Uses pmat;
Function MSort( a: vmatrixptr; col, order: integer ): vmatrixptr;
Function Mexp( ROp: vmatrixptr ): vmatrixptr;
Function Mabs( ROp: vmatrixptr ): vmatrixptr;
Function Mlog( ROp: vmatrixptr ): vmatrixptr;
Function Msqrt( ROp: vmatrixptr ): vmatrixptr;
Function Trace( ROp: vmatrixptr ): double;
Function Helm( n: integer ): vmatrixptr;
Function Index( lower, upper, step: integer ): vmatrixptr;
Function Kron( a, b: vmatrixptr ): vmatrixptr;
Function Det( Datain: vmatrixptr ): double;
Function Cholesky( ROp: vmatrixptr ): vmatrixptr;
Procedure QR( Var ROp, Q, R: vmatrixptr; makeq: boolean );
Function Svd( Var A, U, S, V: vmatrixptr;
makeu, makev: boolean ) : integer;
Function Ginv( ROp: vmatrixptr ): vmatrixptr;
Function Fft( ROp: vmatrixptr; INZEE: boolean ): vmatrixptr;
Function Vec( ROp: vmatrixptr ): vmatrixptr;
Function Diag( ROp: vmatrixptr ): vmatrixptr;
Function Shape( ROp: vmatrixptr; nrows: integer ): vmatrixptr;
Type margins = (ALL,ROWS,COLUMNS);
Function Sum( ROp: vmatrixptr; marg: margins ) : vmatrixptr;
Function Sumsq( ROp: vmatrixptr; marg: margins ): vmatrixptr;
Function Cusum( ROp: vmatrixptr ): vmatrixptr;
Function Mmax( ROp: vmatrixptr; marg: margins ): vmatrixptr;
Function Mmin( ROp: vmatrixptr; marg: margins ): vmatrixptr;
Procedure Crow( Var ROp: vmatrixptr; row: integer; c: double );
Procedure Srow( Var ROp: vmatrixptr; row1, row2: integer );
Procedure Lrow( Var ROp: vmatrixptr; row1, row2: integer; c: double );
Procedure Ccol( Var ROp: vmatrixptr; col: integer; c: double );
Procedure Scol( Var ROp: vmatrixptr; col1, col2: integer );
Procedure Lcol( Var ROp: vmatrixptr; col1, col2: integer; c: double );
(************************************************************)
Implementation
Procedure heapsort( Var a: vmatrixptr );
Var
n, s0, s1, s2, s3, s4, s5, s4m1, s5m1 : integer;
t1, ts : double;
Begin
{ Shell-Metzger sort, PC-World, May 1986 }
n := a^.r;
s0 := n;
s1 := n;
s1 := s1 Div 2;
While s1 > 0 Do Begin
s2 := s0 - s1;
For s3 := 1 To s2 Do Begin
s4 := s3;
While s4 > 0 Do Begin
s5 := s4 + s1;
s4m1 := s4;
s5m1 := s5;
If a^.m( s4m1, 1 ) > a^.m( s5m1, 1 ) Then Begin
t1 := a^.m( s4m1, 1 );
a^.mm( s4m1, 1 )^ := a^.m( s5m1, 1 );
a^.mm( s5m1, 1 )^ := t1;
ts := a^.m( s4m1, 2 );
a^.mm( s4m1, 2 )^ := a^.m( s5m1, 2 );
a^.mm( s5m1, 2 )^ := ts;
s4 := s4 - s1;
End
Else Begin
s4 := 0;
End;
End;
End;
s1 := s1 Div 2;
End;
End; { end heapsort }
Function MSort{(a :vmatrixptr; col, order: integer): vmatrixptr;};
Var
sorted, temp : vmatrixptr;
i,j,t : integer;
Begin
a^.Garbage;
new( sorted, makematrix( a^.r, a^.c ) );
If order <> 0 Then Begin
{ sort column col, then reorder other rows }
new( temp, makematrix( a^.r, 2 ) );
For i := 1 To a^.r Do Begin
temp^.mm( i, 1 )^ := a^.m( i, col );
temp^.mm( i, 2 )^ := i ;
End;
heapsort( temp );
For i := 1 To a^.r Do Begin
For j := 1 To a^.c Do Begin
t := Trunc( temp^.m( i, 2 ) );
sorted^.mm( i, j )^ := a^.m( t, j );
End;
End;
End
Else Begin
{ sort row col, then reorder other columns }
new( temp, makematrix( a^.c, 2 ) );
For i := 1 To a^.c Do Begin
temp^.mm( i, 1 )^ := a^.m( col, i );
temp^.mm( i, 2 )^ := i ;
End;
heapsort( temp );
For i := 1 To a^.c Do Begin
For j := 1 To a^.r Do Begin
t := Trunc( temp^.m( i, 2 ) );
sorted^.mm( j, i )^ := a^.m( j, t );
End;
End;
End;
Dispatch^.Push( sorted );
dispose( temp, killVmatrix );
MSort := Dispatch^.ReturnMat;
End;
Function Mexp{(ROp :vmatrixptr): vmatrixptr;};
Var
temp : vmatrixptr;
i,j : integer;
Begin
{ take exponent of all elements }
ROp^.Garbage;
new( temp, makematrix( ROp^.r, ROp^.c ) );
For i := 1 To ROp^.r Do Begin
For j := 1 To ROp^.c Do
temp^.mm( i, j )^ := exp( ROp^.m( i, j ) );
End;
Dispatch^.Push( temp );
Mexp := Dispatch^.ReturnMat;
End;
Function Mabs{(ROp: vmatrixptr): vmatrixptr;};
Var
temp : vmatrixptr;
i,j : integer;
Begin
{ take absolute values of all elements }
ROp^.Garbage;
new( temp, makematrix( ROp^.r, ROp^.c ) );
For i := 1 To ROp^.r Do Begin
For j := 1 To ROp^.c Do
temp^.mm( i, j )^ := abs( ROp^.m( i, j ) );
End;
Dispatch^.Push( temp );
Mabs := Dispatch^.ReturnMat;
End;
Function Mlog{(ROp :vmatrixptr): vmatrixptr;};
Var
temp : vmatrixptr;
i,j : integer;
R : double;
Begin
{ take logarithm of all elements }
ROp^.Garbage;
new( temp, makematrix( ROp^.r, ROp^.c ) );
For i := 1 To ROp^.r Do Begin
For j := 1 To ROp^.c Do Begin
R := ROp^.m( i, j );
If R <= 0.0 Then
Dispatch^.Nrerror( 'Mlog: zero or negative argument to log' )
Else temp^.mm( i, j )^ := ln( R );
End;
End;
Dispatch^.Push( temp );
Mlog := Dispatch^.ReturnMat;
End;
Function Msqrt{(ROp: vmatrixptr): vmatrixptr;};
Var
temp : vmatrixptr;
i,j : integer;
R : double;
Begin
{ take sqrt of all elements }
ROp^.Garbage;
new( temp, makematrix( ROp^.r, ROp^.c ) );
For i := 1 To ROp^.r Do Begin
For j := 1 To ROp^.c Do Begin
R := ROp^.m( i, j );
If R < 0 Then
Dispatch^.Nrerror( 'Msqrt: zero or negative argument to sqrt' )
Else temp^.mm( i, j )^ := sqrt( R );
End;
End;
Dispatch^.Push( temp );
Msqrt := Dispatch^.ReturnMat;
End;
Function Trace{(ROp: vmatrixptr): double;};
Var
i : integer;
t : double;
Begin
ROp^.Garbage;
t := 0;
If ROp^.r <> ROp^.c Then
Dispatch^.Nrerror( ' Trace: matrix not square in trace' );
For i := 1 To ROp^.r Do t := t + ROp^.m( i, i );
trace := t;
End;
Function Helm {(n: integer): vmatrixptr;};
Var
i, j : integer;
d, den : double;
temp : vmatrixptr;
Begin
new( temp, makematrix( n, n ) );
For i := 1 To n Do Begin
For j := 1 To n Do Begin
d := 1.0 / sqrt( 0.0 + n );
den := d;
If j > 1 Then den := 1.0 / sqrt( 0.0 + j * (j - 1) );
If (j > 1) And (j < i) Then d := 0;
If (j > 1) And (j = i) Then d := - den * ( 0.0 + (j - 1));
If (j > 1) And (j > i) Then d := den;
temp^.mm( i, j )^ := d;
End;
End;
Dispatch^.Push( temp );
helm := Dispatch^.ReturnMat;
End;
Function Index{( lower, upper, step : integer): vmatrixptr;};
Var
cnter, i : integer;
temp : vmatrixptr;
Begin
If step = 0 Then Dispatch^.Nrerror( 'Index: step is zero' );
cnter := 0;
i := lower;
If lower < upper Then Begin
If step < 0 Then
Dispatch^.Nrerror( 'Index: trying to step from lower to upper with negative step size' );
While i <= upper Do Begin
cnter := cnter + 1;
i := i + step;
End;
End
Else Begin
If step > 0 Then
Dispatch^.Nrerror( 'Index: trying to step from upper to lower with positive step size' );
While i >= upper Do Begin
cnter := cnter + 1;
i := i + step;
End;
End;
If cnter = 0 Then
Dispatch^.Nrerror( 'Index: trying to allocate a matrix with zero rows' );
new( temp, makematrix( cnter, 1 ) );
cnter := 1;
i := lower;
If lower < upper Then
While i <= upper Do Begin
temp^.mm( cnter, 1 )^ := i;
i := i + step;
cnter := cnter + 1;
End
Else
While i >= upper Do Begin
temp^.mm( cnter, 1 )^ := i;
i := i + step;
cnter := cnter + 1;
End;
Dispatch^.Push( temp );
Index := Dispatch^.ReturnMat;
End;
Function Kron{(a,b:vmatrixptr): vmatrixptr;};
Var
cc,cr,i,j,k,l : integer;
c : vmatrixptr;
Begin
{ Kroniker's product }
a^.Garbage;
b^.Garbage;
cr := a^.r * b^.r;
cc := a^.c * b^.c;
new( c, makematrix( cr, cc ) );
For i := 1 To a^.r Do Begin
For j := 1 To a^.c Do Begin
For k := 1 To b^.r Do Begin
For l := 1 To b^.c Do Begin
c^.mm( (i - 1) * b^.r + k, (j - 1) * b^.c + l )^ := a^.m( i, j ) * b^.m( k, l );
End;
End;
End;
End;
Dispatch^.Push( c );
Kron := Dispatch^.ReturnMat;
End; (* kron *)
{ private function in unit: used by determinant }
Procedure Pivot( Var Data: vmatrixptr; RefRow: integer;
Var Determ: double; Var DetEqualsZero: boolean );
Var
NewRow, i: integer;
temp : double;
Begin
DetEqualsZero := TRUE;
NewRow := RefRow;
While DetEqualsZero And (NewRow < Data^.r) Do Begin
NewRow := NewRow + 1;
If abs( Data^.m( NewRow, RefRow ) ) > 1.0E - 8 Then Begin
For i := 1 To Data^.r Do Begin
temp := Data^.m( NewRow, i );
Data^.mm( NewRow, i )^ := Data^.m( RefRow, i );
Data^.mm( RefRow, i )^ := temp;
End;
DetEqualsZero := FALSE;
Determ := - Determ; (* row swap adjustment to det *)
End;
End;
End;
Function Det{(Datain : vmatrixptr): double;};
Var
Determ : double;
data : vmatrixptr;
coef : extended;
Row, RefRow, Dimen, i : integer;
DetEqualsZero : boolean;
Begin
Datain^.Garbage;
Determ := 1;
If Datain^.r = Datain^.c Then Begin
Dispatch^.Inclevel;
Dimen := Datain^.r;
new( Data, makematrix( Dimen, Dimen ) );
Data := matequals( Data, Datain );
Row := 0;
RefRow := 0;
DetEqualsZero := FALSE;
While (Not DetEqualsZero) And (RefRow < Dimen - 1) Do Begin
RefRow := RefRow + 1;
If abs( Data^.m( RefRow, RefRow ) ) < 1.0E - 8 Then
Pivot( Data, RefRow, Determ, DetEqualsZero );
If (Not DetEqualsZero) Then
For Row := RefRow + 1 To Dimen Do
If abs( Data^.m( Row, RefRow ) ) > 1.0E - 8 Then Begin
Coef := - Data^.m( Row, RefRow ) / Data^.m( RefRow, RefRow );
For i := 1 To Dimen Do
Data^.mm( Row, i )^ := Data^.m( Row, i ) +
(Coef * Data^.m( RefRow, i ));
End;
Determ := Determ * Data^.m( RefRow, RefRow );
End;
If DetEqualsZero Then Determ := 0
Else Determ := Determ * Data^.m( Dimen, Dimen );
dispose( Data, killvmatrix );
Dispatch^.Declevel;
End
Else Dispatch^.Nrerror( 'Det: Matrix is not square\n' );
det := Determ;
End;
{ / --------------- cholesky decomposition ---------------------- }
{ / ROp is symmetric, returns upper triangular S where ROp := S'S }
{ / ------------------------------------------------------------- }
Function Cholesky{(ROp: vmatrixptr): vmatrixptr;};
Var
nr, i, j, k : integer;
temp : vmatrixptr;
sum : double;
Begin
nr := ROp^.r;
If ROp^.r <> ROp^.c Then
Dispatch^.Nrerror( 'Cholesky: Input matrix not square' );
ROp^.Garbage;
For i := 1 To nr Do Begin
For j := 1 To i - 1 Do
If abs( ROp^.m( i, j ) - ROp^.m( j, i ) ) > 1.0E - 7 Then
Dispatch^.Nrerror( ' Cholesky: Input matrix is not symmetric' );
End;
new( temp, makematrix( nr, nr ) );
temp := matequals( temp, Fill( nr, nr, 0.0 ) );
For i := 1 To nr Do Begin
For j := i To nr Do Begin
sum := 0.0;
For k := 1 To i - 1 Do
sum := sum + temp^.m( k, i ) * temp^.m( k, j );
If i = j Then
If ROp^.m( i, j ) < sum Then
Dispatch^.Nrerror( ' Cholesky: negative pivot' )
Else
temp^.mm( i, j )^ := sqrt( ROp^.m( i, j ) - sum )
Else
If abs( temp^.m( i, i ) ) < 1.0E - 7 Then
Dispatch^.Nrerror( ' Cholesky: zero or negative pivot' )
Else
temp^.mm( i, j )^ := (ROp^.m( i, j ) - sum) / temp^.m( i, i );
End;
End;
Dispatch^.Push( temp );
Cholesky := Dispatch^.ReturnMat;
End;
Procedure QR{(var ROp, Q, R: vmatrixptr; makeq: boolean);};
Var
rr,c,j,i,k: integer;
s, sigma, beta, sum : double;
Begin
ROp^.Garbage;
Q^.Garbage;
R^.Garbage;
Dispatch^.Inclevel;
rr := ROp^.r;
c := ROp^.c;
Q := matequals( Q, ROp );
R := matequals( R, Fill( c, c, 0.0 ) );
For j := 1 To c Do Begin
sigma := 0.0;
For i := j To rr Do
sigma := sigma + Q^.m( i, j ) * Q^.m( i, j );
If abs( sigma ) <= 1.0e - 10 Then
Dispatch^.Nrerror( ' QR: singular X' );
If Q^.m( j, j ) < 0 Then s := - sqrt( sigma )
Else s := sqrt( sigma );
R^.mm( j, j )^ := s;
beta := 1.0 / (s * Q^.m( j, j ) - sigma);
Q^.mm( j, j )^ := Q^.m( j, j ) - s;
For k := j + 1 To c Do Begin
sum := 0.0;
For i := j To rr Do
sum := sum + Q^.m( i, j ) * Q^.m( i, k );
sum := sum * beta;
For i := j To rr Do
Q^.mm( i, k )^ := Q^.m( i, k ) + Q^.m( i, j ) * sum;
R^.mm( j, k )^ := Q^.m( j, k );
End;
End;
If makeq Then Q := matequals( Q, Mult( ROp, Inv( R ) ) );
Dispatch^.Declevel;
End;
{ / ------------------- Singular Valued Decomposition ---------- }
{ / from EISPACK SVD }
{ / ------------------------------------------------------------ }
Function safehypot( a1, b1: extended ): extended;
{ function to get hypotenuse numbers a1 and b1 }
{ where a:=ln(a1) and b:=ln(b1) }
Var
del, sum, a, b : extended;
Begin
a := 2.0 * ln( abs( a1 ) );
b := 2.0 * ln( abs( b1 ) );
If a > b Then del := a
Else del := b;
{ add a1^2 + b1^2, but in logarithms }
sum := ln( exp( a - del ) + exp( b - del ) ) + del;
safehypot := (exp( 0.5 * sum ));
End;
Procedure svdtemp( a: vmatrixptr;
Var w: vmatrixptr;
matu: boolean; Var u: vmatrixptr;
matv: boolean; Var v: vmatrixptr;
Var ierr: integer; Var rv1: vmatrixptr );
Var
m, n, i, j, k, l, ii, i1, kk, k1, ll, l1, its, mn : integer;
c, f, g, h, s, x, y, z, eps, scale, machep, fss : extended;
{
apologies to the purists, but I wanted to preserve the
FORTRAN style here because SVD is a good use of goto's. It
will also help you check it if you go to the original code.
}
Label l190, l210, l270, l290, l360, l390, l410,
l430, l460, l475, l490, l510, l520, l540,
l565, l600, l650, l700, l1000;
Begin
Dispatch^.Inclevel;
a^.Garbage;
w^.Garbage;
u^.Garbage;
v^.Garbage;
rv1^.Garbage;
m := a^.r;
n := a^.c;
{ vmatrixptr a(' a' ,m,n),w(' w' ,n,1), }
{ u(' u' ,m,m),v(' v' ,n,n),rv1(' rv1' ,n,1); }
{ boolean matu,matv; }
{ ********** Machep is a machine dependent parameter specifying }
{ the relative precision of floating point aritmetic. }
{ for pc doubles, range is 1.6^-308 to 1.6^308 }
{ ************** }
machep := exp( - 308.0 * ln( 1.6 ) );
ierr := 0;
u := matequals( u, a );
{ / ********householder reduction to bi diagonal form ******** }
g := 0.0;
scale := 0.0;
x := 0.0;
For i := 1 To n Do Begin
l := i + 1;
rv1^.mm( i, 1 )^ := scale * g;
g := 0.0;
s := 0.0;
scale := 0.0;
If i > m Then Goto l210;
For k := i To m Do scale := scale + abs( u^.m( k, i ) );
If scale = 0.0 Then Goto l210;
For k := i To m Do Begin
u^.mm( k, i )^ := u^.m( k, i ) / scale;
s := s + u^.m( k, i ) * u^.m( k, i );
End;
f := u^.m( i, i );
{g := -((f >= 0.0) ? abs(sqrt(s)) : - abs(sqrt(s)));}
If f >= 0.0 Then g := - abs( sqrt( s ) )
Else g := abs( sqrt( s ) );
h := f * g - s;
u^.mm( i, i )^ := f - g;
If i = n Then Goto l190;
For j := l To n Do Begin
s := 0.0;
For k := i To m Do s := s + u^.m( k, i ) * u^.m( k, j );
f := s / h;
For k := i To m Do u^.mm( k, j )^ := u^.m( k, j ) + f * u^.m( k, i );
End;
l190 :
For k := i To m Do u^.mm( k, i )^ := scale * u^.m( k, i );
l210 : w^.mm( i, 1 )^ := scale * g;
g := 0.0;
s := 0.0;
scale := 0.0;
If (i > m) Or (i = n) Then Goto l290;
For k := l To n Do scale := scale + abs( u^.m( i, k ) );
If scale = 0.0 Then Goto l290;
For k := l To n Do Begin
u^.mm( i, k )^ := u^.m( i, k ) / scale;
s := s + u^.m( i, k ) * u^.m( i, k );
End;
f := u^.m( i, l );
{g := -((f >= 0.0) ? abs(sqrt(s)) : - abs(sqrt(s)));}
If f >= 0.0 Then g := - abs( sqrt( s ) )
Else g := abs( sqrt( s ) );
h := f * g - s;
u^.mm( i, l )^ := f - g;
For k := l To n Do rv1^.mm( k, 1 )^ := u^.m( i, k ) / h;
If i = m Then Goto l270;
For j := l To m Do Begin
s := 0.0;
For k := l To n Do s := s + u^.m( j, k ) * u^.m( i, k );
For k := l To n Do u^.mm( j, k )^ := u^.m( j, k ) + s * rv1^.m( k, 1 );
End;
l270 :
For k := l To n Do u^.mm( i, k )^ := scale * u^.m( i, k );
l290 :
{x := (x > abs(w^.m(i, 1)) + abs(rv1^.m(i, 1))) ? x :
abs(w^.m(i, 1)) + abs(rv1^.m(i, 1)); }
If x <= abs( w^.m( i, 1 ) ) + abs( rv1^.m( i, 1 ) ) Then
x := abs( w^.m( i, 1 ) ) + abs( rv1^.m( i, 1 ) );
End;
{ ********** accumulation of right-hand transformations ******** }
If Not matv Then Goto l410;
For ii := 1 To n Do Begin
i := n + 1 - ii;
If i = n Then Goto l390;
If g = 0.0 Then Goto l360;
{ double division avoids possible underflow }
For j := l To n Do v^.mm( j, i )^ := (u^.m( i, j ) / u^.m( i, l )) / g;
For j := l To n Do Begin
s := 0.0;
For k := l To n Do s := s + u^.m( i, k ) * v^.m( k, j );
For k := l To n Do v^.mm( k, j )^ := v^.m( k, j ) + s * v^.m( k, i );
End;
l360 :
For j := l To n Do Begin
v^.mm( i, j )^ := 0.0;
v^.mm( j, i )^ := 0.0;
End;
l390 : v^.mm( i, i )^ := 1.0;
g := rv1^.m( i, 1 );
l := i;
End;
{*************accumulation of left-hand transformations }
l410 :
If Not matu Then Goto l510;
{ / ************* for i:= min(m,n) step -1 until 1 do ****** }
mn := n;
If m < n Then mn := m;
For ii := 1 To mn Do Begin
i := mn + 1 - ii;
l := i + 1;
g := w^.m( i, 1 );
If i = n Then Goto l430;
For j := l To n Do u^.mm( i, j )^ := 0.0;
l430 :
If g = 0.0 Then Goto l475;
If i = mn Then Goto l460;
For j := l To n Do Begin
s := 0.0;
For k := l To m Do s := s + u^.m( k, i ) * u^.m( k, j );
{ double division avoids possible underflow }
f := (s / u^.m( i, i )) / g;
For k := i To m Do u^.mm( k, j )^ := u^.m( k, j ) + f * u^.m( k, i );
End;
l460 :
For j := i To m Do u^.mm( j, i )^ := u^.m( j, i ) / g;
Goto l490;
l475 :
For j := i To m Do u^.mm( j, i )^ := 0.0;
l490 : u^.mm( i, i )^ := u^.m( i, i ) + 1.0;
End;
{ ********** diagonalization of the bidiagonal form *********** }
l510 : eps := machep * x;
{ ********** for k:=n step -1 until 1 do *********************** }
For kk := 1 To n Do Begin
k1 := n - kk;
k := k1 + 1;
its := 0;
{ ********** test for splitting .********** }
{ for l:=k step -1 until 1 do }
l520 :
For ll := 1 To k Do Begin
l1 := k - ll;
l := l1 + 1;
If abs( rv1^.m( l, 1 ) ) <= eps Then Goto l565;
{ / rv1(1) is always zero, so there is no exit }
{ / through the bottom of the loop ********* }
If abs( w^.m( l1, 1 ) ) <= eps Then Goto l540;
End;
{ ************* cancellation of rv1(l) if l greater than 1****** }
l540 : c := 0.0;
s := 1.0;
For i := l To k Do Begin
f := s * rv1^.m( i, 1 );
rv1^.mm( i, 1 )^ := c * rv1^.m( i, 1 );
If abs( f ) <= eps Then Goto l565;
g := w^.m( i, 1 );
h := safehypot( f, g );
w^.mm( i, 1 )^ := h;
c := g / h;
s := - f / h;
If matu Then Begin
{ go to 560 }
For j := 1 To m Do Begin
y := u^.m( j, l1 );
z := u^.m( j, i );
u^.mm( j, l1 )^ := y * c + z * s;
u^.mm( j, i )^ := - y * s + z * c;
End;
End;
End;
{ ************** test for convergence ******************** }
l565 : z := w^.m( k, 1 );
If l = k Then Goto l650;
{ *************** if shift from bottom 2 by 2 minor ******* }
If its = 50 Then Goto l1000;
its := its + 1;
x := w^.m( l, 1 );
y := w^.m( k1, 1 );
g := rv1^.m( k1, 1 );
h := rv1^.m( k, 1 );
f := ((y - z) * (y + z) + (g - h) * (g + h)) / (2.0 * h * y);
g := safehypot( f, 1.0 );
{fss :=(f >= 0.0) ? abs(g) : - abs(g);}
If f >= 0.0 Then fss := abs( g )
Else fss := - abs( g );
f := ((x - z) * (x + z) + h * (y / (f + fss) - h)) / x;
{ **************** next qr transformation *************** }
c := 1.0;
s := 1.0;
For i1 := l To k1 Do Begin
i := i1 + 1;
g := rv1^.m( i, 1 );
y := w^.m( i, 1 );
h := s * g;
g := c * g;
z := safehypot( f, h );
rv1^.mm( i1, 1 )^ := z;
c := f / z;
s := h / z;
f := x * c + g * s;
g := - x * s + g * c;
h := y * s;
y := y * c;
If matv Then Begin
{ goto l575; }
For j := 1 To n Do Begin
x := v^.m( j, i1 );
z := v^.m( j, i );
v^.mm( j, i1 )^ := x * c + z * s;
v^.mm( j, i )^ := - x * s + z * c;
End;
End;
z := safehypot( f, h );
w^.mm( i1, 1 )^ := z;
{ ************ rotation can be arbitrary if z is zero ******* }
If z <> 0.0 Then Begin
{ goto l580; }
c := f / z;
s := h / z;
End;
{ l580: }
f := c * g + s * y;
x := - s * g + c * y;
If matu Then Begin
{//goto l600; }
For j := 1 To m Do Begin
y := u^.m( j, i1 );
z := u^.m( j, i );
u^.mm( j, i1 )^ := y * c + z * s;
u^.mm( j, i )^ := - y * s + z * c;
End;
End;
l600 : x := x;
End; {//continue}
rv1^.mm( l, 1 )^ := 0.0;
rv1^.mm( k, 1 )^ := f;
w^.mm( k, 1 )^ := x;
Goto l520;
{ / ************** convergence ************** }
l650 :
If z >= 0.0 Then Goto l700;
{ / ************* w(k) is made non-negative ********* }
w^.mm( k, 1 )^ := - z;
If Not matv Then Goto l700;
For j := 1 To n Do v^.mm( j, k )^ := - v^.m( j, k );
l700 : x := x;
End;
ierr := 0;
Dispatch^.Declevel;
exit;
{ ***************** set error -- no convergence to a }
{ singular value after 30 iterations }
l1000 : ierr := k;
Dispatch^.Declevel;
exit;
End;
Function Svd{(var A, U, S, V : vmatrixptr; makeu, makev: boolean;) : integer;};
Var
aa, uu, ss, vv, rv1 : vmatrixptr;
ierr : integer;
Begin
Dispatch^.Inclevel;
A^.Garbage;
ierr := 0;
new( aa, makematrix( 1, 1 ) );
new( uu, makematrix( 1, 1 ) );
new( ss, makematrix( 1, 1 ) );
new( vv, makematrix( 1, 1 ) );
new( rv1, makematrix( 1, 1 ) );
If A^.r < A^.c Then aa := matequals( aa, Tran( A ) )
Else aa := matequals( aa, A );
uu := matequals( uu, Fill( aa^.r, aa^.r, 0.0 ) );
vv := matequals( vv, Fill( aa^.c, aa^.c, 0.0 ) );
ss := matequals( ss, Fill( aa^.c, 1, 0.0 ) );
rv1 := matequals( rv1, Fill( aa^.c, 1, 0.0 ) );
svdtemp( aa, ss, makeu, uu, makev, vv, ierr, rv1 );
If A^.r < A^.c Then Begin
If makeu Then U := matequals( u, vv );
If makev Then V := matequals( v, uu );
End
Else Begin
If makeu Then U := matequals( u, uu );
If makev Then V := matequals( v, vv );
End;
S := matequals( s, ss );
dispose( aa, killvmatrix );
dispose( uu, killvmatrix );
dispose( ss, killvmatrix );
dispose( vv, killvmatrix );
dispose( rv1, killvmatrix );
Dispatch^.Declevel;
svd := ierr;
End;
{ ---------------- end svd ------------------------------------ }
Function Ginv{(ROp : vmatrixptr): vmatrixptr;};
Var
u, s, v, g : vmatrixptr;
i : integer;
t : double;
Begin
ROp^.Garbage;
Dispatch^.Inclevel;
new( u, makematrix( 1, 1 ) );
new( s, makematrix( 1, 1 ) );
new( v, makematrix( 1, 1 ) );
new( g, makematrix( 1, 1 ) );
Svd( ROp, u, s, v, TRUE, TRUE );
For i := 1 To s^.r Do Begin
t := s^.m( i, 1 );
s^.mm( i, 1 )^ := 0;
If abs( t ) <> 0.0 Then s^.mm( i, 1 )^ := 1.0 / t;
End;
g := matequals( g, mult( mult( v, Diag( s ) ), Tran( u ) ) );
dispose( u, killvmatrix );
dispose( s, killvmatrix );
dispose( v, killvmatrix );
Dispatch^.Push( g );
Ginv := Dispatch^.DecReturn;
End;
{ / -------------------------- fft ------------------------------ }
{ / de Boor (1980) SIAM J SCI. STAT. COMPUT. V1 no 1, pp 173-178 }
{ / and NewMat by R. G. Davies a C++ matrix Package }
{ / ------------------------------------------------------------- }
Procedure cossin( n, d: integer; Var c, s: double );
Var
{ calculate cos(twopi*n/d) and sin(twopi*n/d) }
{ minimise roundoff error }
n4 : longint;
sector, sign : integer;
ratio, dn4, dd, divis : double;
Begin
n4 := n * 4;
dn4 := n4;
dd := d;
divis := dn4 / dd;
If divis < 0 Then divis := divis - 0.5
Else divis := divis + 0.5;
sector := trunc( divis );
n4 := n4 - sector * d;
dn4 := n4;
If sector < 0 Then sector := 3 - (3 - sector) Mod 4
Else sector := sector Div 4;
ratio := 1.5707963267948966192 * dn4 / dd;
Case sector Of
0 : Begin
c := cos( ratio ); s := sin( ratio );
End;
1 : Begin
c := - sin( ratio ); s := cos( ratio );
End;
2 : Begin
c := - cos( ratio ); s := - sin( ratio );
End;
3 : Begin
c := sin( ratio ); s := - cos( ratio );
End;
Else writeln( 'got here' );
End;
End;
Procedure fftstep( Var AB, XY: vmatrixptr; after, now, before: integer );
Var
gamma, delta, x, m, j, a, ia, x1, a1, x2, ib, a2, iin : integer;
r_arg, i_arg, r_value, i_value, temp : double;
Begin
gamma := after * before;
delta := now * after;
r_arg := 1.0; i_arg := 0.0;
x := 1;
m := AB^.r - gamma;
For j := 0 To now - 1 Do Begin
a := 1;
x1 := x;
x := x + after;
For ia := 0 To after - 1 Do Begin
{ / generate sins & cosines explicitly rather than iteratively }
{ / for more accuracy; but slower }
cossin( - (j * after + ia), delta, r_arg, i_arg );
a1 := a; a := a + 1;
x2 := x1; x1 := x1 + 1;
If now = 2 Then Begin
ib := before;
While ib <> 0 Do Begin
ib := ib - 1;
a2 := m + a1;
a1 := a1 + after;
r_value := AB^.m( a2, 1 );
i_value := AB^.m( a2, 2 );
XY^.mm( x2, 1 )^ := r_value * r_arg - i_value * i_arg + AB^.m( a2 - gamma, 1 );
XY^.mm( x2, 2 )^ := r_value * i_arg + i_value * r_arg + AB^.m( a2 - gamma, 2 );
x2 := x2 + delta;
End;
End
Else Begin
ib := before;
While ib <> 0 Do Begin
ib := ib - 1;
a2 := m + a1;
a1 := a1 + after;
r_value := AB^.m( a2, 1 );
i_value := AB^.m( a2, 2 );
iin := now - 1;
While iin <> 0 Do Begin
iin := iin - 1;
a2 := a2 - gamma;
temp := r_value;
r_value := r_value * r_arg - i_value * i_arg + AB^.m( a2, 1 );
i_value := temp * i_arg + i_value * r_arg + AB^.m( a2, 2 );
End;
XY^.mm( x2, 1 )^ := r_value;
XY^.mm( x2, 2 )^ := i_value;
x2 := x2 + delta;
End;
End;
End;
End;
End;
Function Fft{(ROp: vmatrixptr; INZEE : boolean): vmatrixptr;};
{ / if INZEE = TTRUE then calculate fft }
{ / if INZEE = FFALSE then calculate inverse fft }
Var
n, nextmx, after, before, next, b1, now, i: integer;
prime : Array[1..26] Of integer;
dn : double;
AB, XY : vmatrixptr;
iinzee : boolean;
Label foundit;
Begin
ROp^.Garbage;
If (ROp^.c <> 1) And (ROp^.c <> 2 ) Then
Dispatch^.Nrerror( 'Fft: Input must have 1 or 2 columns' );
n := ROp^.r; {// length of arrays }
dn := n;
nextmx := 26;
prime[1] := 2; prime[2] := 3; prime[3] := 5; prime[4] := 7;
prime[5] := 11; prime[6] := 13; prime[7] := 17; prime[8] := 19;
prime[9] := 23; prime[10] := 29; prime[11] := 31; prime[12] := 37;
prime[13] := 41; prime[14] := 43; prime[15] := 47; prime[16] := 53;
prime[17] := 59; prime[18] := 61; prime[19] := 67; prime[20] := 71;
prime[21] := 73; prime[22] := 79; prime[23] := 83; prime[24] := 89;
prime[25] := 97; prime[26] := 101;
after := 1;
before := n;
next := 1;
iinzee := TRUE;
Dispatch^.Inclevel;
new( AB, makematrix( n, 2 ) );
new( XY, makematrix( n, 2 ) );
If ROp^.c = 1 Then Begin
For i := 1 To n Do Begin
AB^.mm( i, 1 )^ := ROp^.m( i, 1 );
AB^.mm( i, 2 )^ := 0.0;
End;
End
Else AB := matequals( AB, ROp );
XY := matequals( XY, Fill( n, 2, 0.0 ) );
If Not INZEE Then Begin
{ take complex conjugate for ifft }
For i := 1 To n Do AB^.mm( i, 2 )^ := - AB^.m( i, 2 );
End;
Repeat
While TRUE Do Begin
If next <= nextmx Then now := prime[next];
b1 := before Div now;
If b1 * now = before Then Goto foundit;
next := next + 1;
now := now + 2;
End;
foundit : before := b1;
If iinzee = TRUE Then fftstep( AB, XY, after, now, before )
Else fftstep( XY, AB, after, now, before );
{iinzee :=(iinzee = TTRUE) ? FFALSE : TTRUE;}
If iinzee = TRUE Then iinzee := FALSE
Else iinzee := TRUE;
after := after * now;
Until before = 1 ;
If iinzee = TRUE Then XY := matequals( XY, AB );
{ divide by n for ifft }
If Not INZEE Then
For i := 1 To XY^.r Do Begin
XY^.mm( i, 1 )^ := XY^.m( i, 1 ) / dn;
XY^.mm( i, 2 )^ := XY^.m( i, 2 ) / dn;
End;
Dispatch^.Push( XY );
dispose( AB, killvmatrix );
fft := Dispatch^.DecReturn;
End;
{ /////////////////// reshaping functions }
Function Vec{( ROp : vmatrixptr): vmatrixptr;};
Var
lrc, lr, lc : longint;
r,c,cnter, i,j : integer;
temp : vmatrixptr;
Begin
{// send columns of ROp to a vector}
ROp^.Garbage;
lr := r;
lc := c;
lrc := ROp^.r * ROp^.c;
If (lrc >= 32768) Or (lrc < 1) Then
Dispatch^.Nrerror( ' Vec: vec produces invalid indices' );
r := ROp^.r * ROp^.c;
c := 1;
new( temp, makematrix( r, c ) );
cnter := 1;
For j := 1 To ROp^.c Do
For i := 1 To ROp^.r Do Begin
temp^.mm( cnter, 1 )^ := ROp^.m( i, j );
cnter := cnter + 1;
End;
Dispatch^.Push( temp );
vec := Dispatch^.ReturnMat;
End;
Function Diag{(ROp : vmatrixptr) : vmatrixptr;};
Var
temp : vmatrixptr;
i, rr, cc : integer;
Begin
{ make a diagonal matrix from a vector or the principal diagonal of}
{ a square matrix, off diagonal elements are zero }
ROp^.Garbage;
If (ROp^.r <> ROp^.c) And (ROp^.c <> 1) Then
Dispatch^.Nrerror( 'Diag: ROp is not square or is not a column vector' );
Dispatch^.Inclevel;
new( temp, makematrix( 1, 1 ) );
temp := matequals( temp, Fill( ROp^.r, ROp^.r, 0.0 ) );
rr := ROp^.r;
cc := ROp^.c;
For i := 1 To rr Do
If rr = cc Then temp^.mm( i, i )^ := ROp^.m( i, i )
Else temp^.mm( i, i )^ := ROp^.m( i, 1 );
Dispatch^.Push( temp );
Diag := Dispatch^.DecReturn;
End;
Function Shape{(ROp: vmatrixptr; nrows: integer) : vmatrixptr;};
Var
nr, lr, lc, lrc : longint;
cnter, i,j : integer;
temp1, temp: vmatrixptr;
Begin
{reshape a matrix into n rows, nrows must divide r*c without a}
{remainder }
ROp^.Garbage;
nr := nrows;
lr := ROp^.r;
lc := ROp^.c;
If (lr * lc Mod nr) <> 0 Then
Dispatch^.Nrerror( ' Shape: nrows divides r*c with a remainder' );
lrc := (lr * lc) Div nr;
If (nr >= 32768) Or (nr < 1) Or (lrc >= 32768) Or (lrc < 1) Then
Dispatch^.Nrerror( ' Shape: reshape produces invalid indices' );
Dispatch^.Inclevel;
new( temp1, makematrix( 1, 1 ) );
temp1 := matequals( temp1, Vec( ROp ) );
new( temp, makematrix( nrows, lrc ) );
cnter := 1;
For i := 1 To temp^.r Do
For j := 1 To temp^.c Do Begin
temp^.mm( i, j )^ := temp1^.m( cnter, 1 );
cnter := cnter + 1;
End;
Dispatch^.Push( temp );
dispose( temp1, killvmatrix );
Shape := Dispatch^.DecReturn;
End;
Function Sum{( ROp : vmatrixptr; marg : margins ) : vmatrixptr;};
Var
i,j,r,c : integer;
ssum : double;
temp : vmatrixptr;
Begin
{ sum(ROp, ALL) returns 1x1}
{ sum(ROp, ROWS) returns 1xc }
{ sum(ROp, COLUMNS) returns rx1 }
ROP^.Garbage;
Case marg Of
ALL : Begin
r := 1; c := 1;
End;
ROWS : Begin
r := 1; c := ROP^.c;
End;
COLUMNS : Begin
r := ROP^.r; c := 1;
End;
Else
Dispatch^.Nrerror( ' Sum: invalid margin specification' );
End;
new( temp, makematrix( r, c ) );
Case marg Of
ALL : Begin
ssum := 0.0;
For i := 1 To ROp^.r Do
For j := 1 To ROp^.c Do ssum := ssum + ROp^.m( i, j );
temp^.mm( 1, 1 )^ := ssum;
End;
ROWS : Begin
For j := 1 To ROp^.c Do Begin
ssum := 0.0;
For i := 1 To ROp^.r Do ssum := ssum + ROp^.m( i, j );
temp^.mm( 1, j )^ := ssum;
End;
End;
COLUMNS : Begin
For i := 1 To ROp^.r Do Begin
ssum := 0.0;
For j := 1 To ROp^.c Do ssum := ssum + ROp^.m( i, j );
temp^.mm( i, 1 )^ := ssum;
End;
End;
End;
Dispatch^.Push( temp );
Sum := Dispatch^.ReturnMat;
End;
Function Sumsq{(ROp: vmatrixptr; marg: margins): vmatrixptr;};
Var
i,j,r,c : integer;
sum : double;
temp : vmatrixptr;
Begin
{ sumsq(ROp, ALL) returns 1x1}
{ sumsq(ROp, ROWS) returns 1xc }
{ sumsq(ROp, COLUMNS) returns cx1 }
ROp^.Garbage;
Case marg Of
ALL : Begin
r := 1; c := 1;
End;
ROWS : Begin
r := 1; c := ROp^.c;
End;
COLUMNS : Begin
r := ROp^.r; c := 1;
End;
Else
Dispatch^.Nrerror( ' Sumsq: invalid margin specification' );
End;
new( temp, makematrix( r, c ) );
Case marg Of
ALL : Begin
sum := 0.0;
For i := 1 To ROp^.r Do
For j := 1 To ROp^.c Do sum := sum + sqr( ROp^.m( i, j ) );
temp^.mm( 1, 1 )^ := sum;
End;
ROWS : Begin
For j := 1 To ROp^.c Do Begin
sum := 0.0;
For i := 1 To ROp^.r Do
sum := sum + sqr( ROp^.m( i, j ) );
temp^.mm( 1, j )^ := sum;
End;
End;
COLUMNS : Begin
For i := 1 To ROp^.r Do Begin
sum := 0.0;
For j := 1 To ROp^.c Do sum := sum + sqr( ROp^.m( i, j ) );
temp^.mm( i, 1 )^ := sum;
End;
End;
End;
Dispatch^.Push( temp );
sumsq := Dispatch^.ReturnMat;
End;
Function Cusum{( ROp: vmatrixptr): vmatrixptr;};
Var
i,j : integer;
sum : double;
temp : vmatrixptr;
Begin
{// cumulative sum along rows}
ROp^.Garbage;
new( temp, makematrix( ROp^.r, ROp^.c ) );
sum := 0.0;
For i := 1 To ROp^.r Do Begin
For j := 1 To ROp^.c Do Begin
sum := sum + ROp^.m( i, j );
temp^.mm( i, j )^ := sum;
End;
End;
Dispatch^.Push( temp );
Cusum := Dispatch^.ReturnMat;
End;
Function Mmax{( ROp: vmatrixptr; marg : margins): vmatrixptr;};
Var
i,j,r,c,maxr,maxc : integer;
mmmax : double;
temp : vmatrixptr;
Begin
{// Mmax(ROp, ALL) returns 3x1}
{/ Mmax(ROp, ROWS) returns 3xc }
{/ Mmax(ROp, COLUMNS) returns rx3 }
ROp^.Garbage;
Case marg Of
ALL : Begin
r := 3; c := 1;
End;
ROWS : Begin
r := 3; c := ROp^.c;
End;
COLUMNS : Begin
r := ROp^.r; c := 3;
End;
Else
Dispatch^.Nrerror( ' Mmax: invalid margin specification' );
End;
new( temp, makematrix( r, c ) );
Case marg Of
ALL : Begin
mmmax := ROp^.m( 1, 1 );
maxr := 1;
maxc := 1;
For i := 1 To ROp^.r Do
For j := 1 To ROp^.c Do
If mmmax < ROp^.m( i, j ) Then Begin
mmmax := ROp^.m( i, j );
maxr := i;
maxc := j;
End;
temp^.mm( 1, 1 )^ := maxr;
temp^.mm( 2, 1 )^ := maxc;
temp^.mm( 3, 1 )^ := mmmax;
End;
ROWS : Begin
For j := 1 To c Do Begin
mmmax := ROp^.m( 1, j );
maxr := 1;
maxc := j;
For i := 1 To ROp^.r Do
If mmmax < ROp^.m( i, j ) Then Begin
maxr := i;
mmmax := ROp^.m( i, j );
End;
temp^.mm( 1, j )^ := maxr;
temp^.mm( 2, j )^ := maxc;
temp^.mm( 3, j )^ := mmmax;
End;
End;
COLUMNS : Begin
For i := 1 To r Do Begin
mmmax := ROp^.m( i, 1 );
maxr := i;
maxc := 1;
For j := 1 To ROp^.c Do
If mmmax < ROp^.m( i, j ) Then Begin
maxc := j;
mmmax := ROp^.m( i, j );
End;
temp^.mm( i, 1 )^ := maxr;
temp^.mm( i, 2 )^ := maxc;
temp^.mm( i, 3 )^ := mmmax;
End;
End;
End;
Dispatch^.Push( temp );
Mmax := Dispatch^.ReturnMat;
End;
Function Mmin{( ROp: vmatrixptr; marg : margins): vmatrixptr;};
Var
i,j,r,c,minr,minc : integer;
mmmin : double;
temp : vmatrixptr;
Begin
{/ Mmin(ROp, ALL) returns 3x1}
{/ Mmin(ROp, ROWS) returns 3xc }
{/ Mmin(ROp, COLUMNS) returns rx3 }
ROp^.Garbage;
Case marg Of
ALL : Begin
r := 3; c := 1;
End;
ROWS : Begin
r := 3; c := ROp^.c;
End;
COLUMNS : Begin
r := ROp^.r; c := 3;
End;
Else
Dispatch^.Nrerror( ' Mmin: invalid margin specification' );
End;
new( temp, makematrix( r, c ) );
Case marg Of
ALL : Begin
mmmin := ROp^.m( 1, 1 );
minr := 1;
minc := 1;
For i := 1 To ROp^.r Do
For j := 1 To ROp^.c Do
If mmmin > ROp^.m( i, j ) Then Begin
mmmin := ROp^.m( i, j );
minr := i;
minc := j;
End;
temp^.mm( 1, 1 )^ := minr;
temp^.mm( 2, 1 )^ := minc;
temp^.mm( 3, 1 )^ := mmmin;
End;
ROWS : Begin
For j := 1 To c Do Begin
mmmin := ROp^.m( 1, j );
minr := 1;
minc := j;
For i := 1 To ROp^.r Do
If mmmin > ROp^.m( i, j ) Then Begin
minr := i;
mmmin := ROp^.m( i, j );
End;
temp^.mm( 1, j )^ := minr;
temp^.mm( 2, j )^ := minc;
temp^.mm( 3, j )^ := mmmin;
End;
End;
COLUMNS : Begin
For i := 1 To r Do Begin
mmmin := ROp^.m( i, 1 );
minr := i;
minc := 1;
For j := 1 To ROp^.c Do
If mmmin > ROp^.m( i, j ) Then Begin
minc := j;
mmmin := ROp^.m( i, j );
End;
temp^.mm( i, 1 )^ := minr;
temp^.mm( i, 2 )^ := minc;
temp^.mm( i, 3 )^ := mmmin;
End;
End;
End;
Dispatch^.Push( temp );
Mmin := Dispatch^.ReturnMat;
End;
{ /////////////////////////////// elemenatary row and column operations }
{ ///////////////// note these functions operate directly on ROp }
Procedure Crow{(var ROp: vmatrixptr; row : integer; c : double)};
Var
i : integer;
Begin
{// multiply a row by a constant}
ROp^.Garbage;
If (row < 1) Or (row > ROp^.r) Then
Dispatch^.Nrerror( ' Crow: row index out of range' );
For i := 1 To ROp^.c Do ROp^.mm( row, i )^ := c * ROp^.m( row, i );
End;
Procedure Srow{(var ROp : vmatrixptr; row1, row2 : integer);};
Var
i,j : integer;
tmp : double;
Begin
{// swap rows }
ROp^.Garbage;
If (row1 < 1) Or (row1 > ROp^.r) Or (row2 < 1) Or (row2 > ROp^.r) Then Dispatch^.Nrerror( ' Srow: row index out of range' );
If row1 = row2 Then exit;
For i := 1 To ROp^.c Do Begin
tmp := ROp^.m( row1, i );
ROp^.mm( row1, i )^ := ROp^.m( row2, i );
ROp^.mm( row2, i )^ := tmp;
End;
End;
Procedure Lrow{(var ROp : vmatrixptr; row1, row2 : integer; c : double);};
Var
i : integer;
Begin
{// add c*r1 to r2}
ROp^.Garbage;
If (row1 < 1) Or (row1 > ROp^.r) Or (row2 < 1) Or (row2 > ROp^.r) Then Dispatch^.Nrerror( ' Lrow: row index out of range' );
For i := 1 To ROp^.c Do
ROp^.mm( row2, i )^ := ROp^.m( row2, i ) + c * ROp^.m( row1, i );
End;
Procedure Ccol{(var ROp : vmatrixptr; col : integer; c : double);};
Var
i : integer;
Begin
{// multiply a col by a constant}
ROp^.Garbage;
If (col < 1) Or (col > ROp^.c) Then
Dispatch^.Nrerror( ' Ccol: col index out of range' );
For i := 1 To ROp^.r Do ROp^.mm( i, col )^ := c * ROp^.m( i, col );
End;
Procedure Scol{(var ROp : vmatrixptr; col1, col2 : integer);};
Var
i : integer;
tmp : double;
Begin
{// swap cols}
ROp^.Garbage;
If (col1 < 1) Or (col1 > ROp^.c) Or (col2 < 1) Or (col2 > ROp^.c) Then Dispatch^.Nrerror( ' Scol: col index out of range' );
If col1 = col2 Then exit;
For i := 1 To ROp^.r Do Begin
tmp := ROp^.m( i, col1 );
ROp^.mm( i, col1 )^ := ROp^.m( i, col2 );
ROp^.mm( i, col2 )^ := tmp;
End;
End;
Procedure Lcol{(var ROp : vmatrixptr; col1, col2 : integer; c : double);};
Var
i: integer;
Begin
{// add c*c1 to c2}
ROp^.Garbage;
If (col1 < 1) Or (col1 > ROp^.c) Or (col2 < 1) Or (col2 > ROp^.c) Then Dispatch^.Nrerror( ' Lcol: col index out of range' );
For i := 1 To ROp^.r Do
ROp^.mm( i, col2 )^ := ROp^.m( i, col2 ) + c * ROp^.m( i, col1 );
End;
Begin
End.